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£\j , A method is proposed to compute the interfacial free energy of a Lennard- Jones system in contact 

with a structured wall by molecular dynamics simulation. Both the bulk liquid and bulk face- 
centered-cubic crystal phase along the (111) orientation are considered. Our approach is based 
£Nj , on a thermodynamic integration scheme where first the bulk Lennard- Jones system is reversibly 

» ( ■ transformed to a state where it interacts with a structureless flat wall. In a second step, the flat 

f"S |' structureless wall is reversibly transformed into an atomistic wall with crystalline structure. The 

dependence of the interfacial free energy on various parameters such as the wall potential, the 
density and orientation of the wall is investigated. The conditions are indicated under which a 
\Q ' Lennard- Jones crystal partially wets a fiat wall. 
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I. INTRODUCTION 

Knowledge of the interfacial free energy between a 
crystal or liquid in contact with a solid wall is crucial 
to the understanding of heterogeneous nucleation and 
wetting phenomena [lH5(. However, interfacial free en- 
ergies are hardly accessible in experiments and in fact 
only a few measurements have been reported so far (see 
e.g. [111,0). 

Due to the lack of experimental data, particle-based 
simulation techniques such as Molecular Dynamics (MD) 
and Monte Carlo (MC) are of special importance to 
understand the properties of wall-liquid and wall-crystal 
interfaces and to rationalize calculations in the frame- 
work of density functional theory (l(il4l2l |. In this con- 
text, MC and MD simulations have been used to un- 
derstand the microscopic mechanism of fluid wetting on 
solid surfaces d [l3l - [Tq | as well as the wetting and dry- 
ing transition of a fluid at liquid- vapor coexistence and 
in contact with a solid wall [15l - [l8j . The question of how 
the wall structure affects the interfacial tension with re- 
spect to liquid , vapor and solid phases has been also ad- 
dressed [l6l fl9|| . 

On a macroscopic scale, a crystal that partially wets 
a wall might be described as a spherical cap. Then, the 
contact angle 6 C of the cap with the wall is given by 
Young's equation [5j, 



7wc + Yd cos 6 C = 7wi 



(1) 



with 7 WC the wall-crystal, y c \ the crystal-liquid, and 7 w i 
the wall-liquid interfacial free energy. Equation ([1]) de- 
scribes the condition of a spherical crystal droplet resting 
on a wall, being in coexistence with the liquid phase. In- 
complete wetting corresponds to contact angles < 6 C < 

IT. 

On a nanoscopic scale, deviations from Young's equa- 
tion can be expected, e.g . due to the contribution of line 
tension effects la, [20|. To quantify the latter devia- 
tions, reliable estimates of 7 WC , y c \ and 7 w i are required. 
Then, the contact angle can be obtained via Eq. ([T]) and 



compared to a direct measurement of 9 C . 

In this paper, we propose a thermodynamic integration 
(TI) scheme for the calculation of 7 WC and 7 w j. To 
obtain 7 w j, most previous studies have used the mechan- 
ical approach of calculating the normal and tangential 
pressure components at the wall and integrating over the 
pressure anisotropy (PA) 0, El EMI . While the PA 
method is valid for planar wall-liquid or liquid- vapor in- 
terfaces it fails in case of small liquid drops in contact 
with a solid wall [HJ. Moreover, its use is justified only 
for systems where the interfacial tension equals the inter- 
facial free energy [23[ . This is true for a wall represented 
by a time- independent external field [H HU, HH or a wall 
made of particles rigidly fixed at the sites of an ideal 
lattice [l5rll8l [H [27|. However, for systems which can 
support stress, such as a wall consisting of a "fully in- 
teracting solid phase" [28| . this method is invalid. For 
the same reason, the PA technique cannot be used to de- 
termine 7 W c [HI]. Even for wall-liquid interfaces, the PA 
method can yield results with acceptable precision only 
with huge computational effort. Most previous works 
based on the PA technique yielded results of low accu- 
racy and the values of the interfacial tension reported in 
the literature differ widely, even for simple systems. 

Due to the obvious disadvantage in using the PA 
method, a few thermodynamic approaches have been de- 
veloped to evaluate the wall-liquid and wall-crystal in- 
terfacial free energies with improved precision. Heni and 
Lowen [2t| combined MC simulations and thermody- 
namic integration to determine the interfacial free en- 
ergies of hard sphere liquids and solids near a planar 
structureless wall over a whole range of bulk densities 
including the solid-liquid coexistence density. In their 
thermodynamic integration scheme, a bulk hard sphere 
system was reversibly transformed into a system interact- 
ing with a more and more impenetrable wall and finally 
a hard wall. Fortini and Djikstra [3(| used a thermo- 
dynamic integration scheme based on exponential poten- 
tials to calculate 7 w i and 7 WC at bulk coexistence condi- 
tions. Their results were in good agreement with those 
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of Hcni and Lowen but obtained with significantly higher 
precision. Due to precrystallization of the hard spheres 
near the wall close to the bulk freezing transition, both 
Heni and Lowen and Fortini and Dijkstra extrapolated 
the value of the interfacial tensions at coexistence from 
the data at lower densities. 



Laird and Davidchack 31| developed a TI method by 
the use of "cleaving potentials" , to obtain 7 w i and 7 W r, for 
hard sphere systems at coexistence. In another work [28j , 
they used the "Gibbs-Cahn integration" method, to ob- 
tain wall-fluid interfacial free energies for hard sphere 
systems. This method yielded results consistent with 
the TI method with "cleaving potentials" but were ob- 
tained with significantly less computational effort. How- 
ever, "Gibbs-Cahn integration" requires that one knows 
already the interfacial free energy at one point. Deb et 
al. [32|, [33| compared different methods to obtain wall- 
fluid and wall-crystal interfacial free energies for hard 
sphere systems confined by hard walls or, soft walls de- 
scribed by the Weeks- Chandler- Anderson (WCA) poten- 
tial. They introduced a scheme similar to Wang-Landau 
sampling [34| . known as the "ensemble mixing" method, 
to perform a TI from a system without walls to a sys- 
tem confined by walls. For hard spheres, Deb et al. 
obtained good agreement with the results of Laird and 
Davidchack. 

In contrast to these few works on hard-sphere systems, 
there is a dearth of results on the interfacial free energies 
of systems with continuous potentials, such as Lennard- 
Jones (LJ) systems. Recently, Leroy et al. [lj| obtained 
7wi for a LJ liquid in contact with a flexible LJ struc- 
tured wall by the use of a TI technique, known as the 
"phantom wall" method. In this approach, the struc- 
tured wall interacting with the liquid is gradually moved 
away from the liquid, while a structureless flat wall is 
moved towards it such that in the final state, the liquid 
interacts only with the structureless wall. Computing the 
free energy difference during this transformation, along 
with the interfacial free energy of liquid in contact with 
the structureless wall, gives 7 w j. 7 w i for the liquid-flat 
wall system, which serves as the reference state for their 
system was obtained using the PA technique. Since the 
PA technique fails in case of crystal-wall interfaces [23[ 
one cannot use their scheme to determine 7 WC for crystal 
in contact with a structured wall. In fact, much less is 
known about 7 WC for LJ systems in contact with a wall. 

Grochola et al. [3|| developed another TI technique 
which they have called "A-integration" , to determine the 
surface free energies of solids. In principle, this technique 
could be also applied to wall-crystal or wall-liquid inter- 
faces, but the method has not been worked out yet for 
such interfaces. 

A straightforward and comprehensive method is thus 
needed to compute the interfacial free energies of L J sys- 
tems in contact with a wall. In the present work, a novel 
TI scheme is introduced to compute the interfacial free 
energy of a LJ system confined between walls. We con- 
sider both the liquid as well as the fee crystal phase along 



the (HI) orientation near the wall. While most previous 
works employing TI methods to obtain 7 w i or 7 WC are 
limited to structureless walls, here we specifically con- 
sider the case of a structured wall, consisting of particles 
rigidly attached to the sites of an ideal fee lattice. Our 
scheme consists of TI in two steps, providing a reversible 
thermodynamic path that transforms the bulk LJ system 
into a LJ fluid or crystal interacting with a structured 
wall. In the first step, a thermodynamic path is devised 
to reversibly transform the bulk LJ system without walls 
and periodic boundary conditions in all directions to a 
state where it interacts with the structureless wall. This 
is accomplished by gradually modifying the interaction 
potential between the wall and the LJ particles along the 
thermodynamic path. The technique is inspired by the 
method proposed by Heni and Lowen [29j to compute the 
interfacial free energy of hard sphere fluids and crystals 
in contact with a hard wall. 

The LJ system interacting with the flat wall serves as 
the reference state to calculate the interfacial free energy 
of the LJ liquid or crystal in contact with a structured 
wall. In the second step, another TI scheme reversibly 
changes the structureless wall interacting with the LJ 
system into a structured wall. This is done by gradually 
switching off the flat walls and simultaneously switching 
on the structured walls. While previous methods based 
on TI techniques for the calculation of 7 w i make use of 
"cleaving potentials" [3lJ or "phantom walls " [l9[ , here 
we directly modify the interaction potential to make the 
transformation from the reference state to the final state 
in each of the two steps. Though this TI scheme is specif- 
ically developed for a LJ potential, it can be easily gen- 
eralized to more complex potentials. 

The wetting behavior of a liquid or crystal in contact 
with a structured solid wall will be affected by various 
parameters. In this work we focus on three parameters: 
i) the interaction strength between the wall and the LJ 
system, ii) the density of the structured wall, and iii) 
the orientation of the structured wall with respect to the 
interface normal. For these cases, wall-crystal interfacial 
free energies only for the (111) orientation of the crystal 
are considered. 

Since the PA method has been widely applied in the 
past to evaluate the wall-liquid interfacial free energy, we 
will compare results obtained from it with those yielded 
by the TI method for both flat and structured walls. In 
addition, we will also show that the interfacial free ener- 
gies of the LJ system interacting with a flat wall can be 
obtained directly at coexistence, without any extrapola- 
tion from data at low densities, enabling us to investigate 
its wetting behavior. 

Furthermore, to compare the estimates of interfacial 
free energy yielded by our TI scheme with that obtained 
in previous works, we apply our technique to a model 
system studied by Tang and Harris (TH) [16[ using the 
mechanical definition of the interfacial free energy. Their 
system consisted of a LJ fluid confined between identical 
rigid structured walls oriented along the (100) orienta- 
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tion, under conditions of liquid- vapor coexistence. Later, 
Grzelak and Errington (GE) [36jj investigated the same 
system using Grand Canonical Transition Matrix Monte 
Carlo (GCTMMC) simulations. They computed the in- 
terracial free energy profile as a function of the surface 
density at bulk liquid-vapor saturation condition, to ob- 
tain the contact angle and the solid- vapor and solid-liquid 
interfacial tensions. For this system, we will examine the 
variation of 7 w i as a function of the wall-liquid interac- 
tion strength and compare estimates of 7 w j from these 
two studies. Due to the paucity of studies on the crystal- 
wall interfacial free energy, we will restrict this compari- 
son with previous works only to the wall-liquid interfacial 
free energy. 

In the following, we introduce the details of the model 
potentials considered in this work (Sec. [IT]), give the vari- 
ous definitions of interfacial free energies, outline the PA 
method, describe the proposed TI scheme, and provide 
the main details of the simulation (Sec. IIII|) . Then, we 
present the results (Sec. IIVI) and finally draw some con- 
clusions fSec. IVf. 



II. MODEL POTENTIAL 



z direction the particles are confined by the structured 
wall, between z = z^ at the top and z — z t at the bottom. 
The system thus consists of two planar wall-liquid (or 
wall-crystal) interfaces with a total area of A = 2L x L y . 
The structured wall is arranged in a manner such that the 
wall layers closest to the LJ system are positioned at z^ = 
— L z /2 and zt = L z /2. Also, an integer number of unit 
cells was chosen for the structured wall such that the wall 
is exactly adapted to the lateral size of the simulation cell. 

The width of the structured wall is chosen large enough 
to avoid LJ particles on opposite sides of the wall from 
interacting with each other since the determination of 
interfacial free energy by TI or PA methods is built on 
the assumption of two independent wall-liquid (or wall- 
crystal) interfaces. 

The TI scheme adopted in this work consists of two 
steps. First, a bulk LJ system with periodic boundary 
conditions is transformed into a state where the LJ sys- 
tem interacts with impenetrable flat walls. Then, in the 
second step, the flat walls are reversibly transformed into 
structured walls. The structureless flat wall (fw) is taken 
to be a purely repulsive potential interacting along the 
z direction with the LJ particles and is described by a 
WCA potential, 



The MD system to determine the interfacial excess free 
energy of a LJ system in contact with a structured wall 
consists of N identical particles interacting with each 
other and with the structured wall via a shifted-force 
LJ Q potential. If two particles i and j of types a and /3 
are separated by a distance ry , the interaction potential 
is written as 



^ a p(nj) - (j) a p(r c ) - <j>' a p{rij = r c )[ riJ - r c 







for < Tij < r c 
for > r c , 



(2) 

where the prime in <fi' a p denotes the derivative with re- 
spect to r and 
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(3) 



Ufw{z) = { for < z < z cw , 
.0 for z > z cw 



x w(z) 



(4) 



with the cut-off z cw = 2 ' a pp and z = Zi — Z the distance 
of particle i at Zi to one of the flat walls at Z = z\> or 
Z = z t . The function w{z) ensures that Uf w (z) goes 
smoothly to zero at z — z cw and is given by 



w(z) 



1 



h± + (i/(z- Zcw yy 



(5) 



where the dimensionless parameter h is set to 0.005. 

To compare to the results of Tang and Harris |16| and 
Grzelak and Errington [361 ] , we also consider a truncated 
and shifted L J potential for the particle-particle (pp) and 
particle-structured wall (pw) interactions, 



In Eq. @, a or f3 can represent a LJ particle (p) or a 
structured wall particle (w). The parameters e a p and 
a a p have units of energy and length, respectively. The 
cut-off distance is set to r c — 2.5cr Q , ) g. 

In the following, energies, lengths and masses are given 
in units of e pp , a pp and m p , respectively. Thus, tempera- 
ture, pressure and interfacial free energy are expressed in 
units of Cpp/ks, epp/c pp and e pp /tr pp , respectively. Time 
is made dimensionless by reducing it with respect to the 

characteristic time scale J m p <7 2 p /e pp . For simplicity, we 
choose cr wp = er pp . 

The N identical liquid or crystal particles are enclosed 
within a simulation box of size L x xL y xL z , using periodic 
boundary conditions in the x and y directions. In the 



* / x J <t><xp{rn) ~ <t> a p{r c ) for r tJ < r c , 
M ^ (rij) -\0 for r«>r c (6) 

with a/3 = pp, pw and the cut-off radius r € — 2.5cr pp . 
Moreover, we choose er pw = l.lcr pp and vary the parame- 
ter e pw in units of e pp in order to determine the interfacial 
free energy 7 wl as a function of the strength of the pw 
interactions. As Tang and Harris QTij . we use a substrate 
consisting of three layers of atoms rigidly fixed to fee lat- 
tice sites, with the (100) orientation of the wall facing the 
liquid along the z direction. The average number density 
of the liquid is set to p p — 0.661 cr pp 3 and that of the 
substrate to p w = 0.59 cr pp keeping the temperature of 
the system fixed at T = 0.9/cb/£ pp . 
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III. CALCULATION OF INTERFACIAL FREE 
ENERGIES 

A. Definitions 

The Hamiltonian of our model, corresponding to the 
LJ system interacting with a solid wall, can be written 

as 

N p 1 W p N p 

H(r, P) = J2 +EE m ppM + ^waii (7) 



i=i 



i=l j=i+l 



with N p the total number of LJ particles and U wa n the 
wall-particle potential. For interactions of the LJ system 
with a flat wall, U wa \\ = X)i=i u fv( z = z i ~ Z), and for 
the system in contact with a structured wall, [Avail = 
Si^i Sj=i u pw(r ij) (with iV w the total number of wall 
particles). In Eq. ([7]). there is no kinetic energy term 
for the walls, since the fiat walls are considered to be 
of infinite mass and immovable; similarly, the structured 
wall particles are considered to be immobile. 

Our simulations are performed in the NP^AT ensem- 
ble, where the number of particles N, surface area A and 
temperature T are kept constant and the length of the 
simulation box along the z direction is allowed to fluc- 
tuate in order to maintain a constant normal pressure 
Pn- The use of the NPjyAT ensemble is necessary to 
maintain a constant bulk density of the system when TI 
is applied (see below). Moreover, any stress present in 
the crystal due to interaction with the walls can relax 
during the NP^AT simulation. The determination of 
the interfacial free energy by thermodynamic or mechan- 
ical approaches demands that there is a bulk region in the 
middle of the simulation box where the density is equal to 
the bulk density of the homogeneous system. Hence, the 
system size along the z direction must be large enough 
to prevent the two walls on cither side of the LJ system 
from influencing each other. 

The isothermal-isobaric partition function correspond- 
ing to the Hamiltonian ([7]) is 



f 



!NP N AT 



h 3N Nl _ 

x AdL z dr N dp N 



exp 



tf(r,p) + P N AL z 



(8) 



where r and p denote respectively the positions and 
momenta of the particles and h is the Planck con- 
stant. The Gibbs free energy G of the confined liq- 
uid or crystal is related to the partition function ((SJ by 

G = —kBTlllQNp N AT- 

The derivative of Gibbs free energy with respect to the 
surface area defines the interfacial tension: 



7 = 



dG 
dA 



(9) 



This thermodynamic definition of the interfacial tension 
is equivalent to the mechanical definition [37| : 



f 



= 2 I PnW-^W 



(10) 



where Pn(z) and Pt{z) are respectively the normal and 
tangential pressure profiles of the liquid and the factor 
1/2 is introduced to account for the fact that the liquid 
is confined between two identical walls. The local pres- 
sure tensor components Pn(z) and Pt(z) are defined in 
Eqs. ([13]). (|T6f and (14]) (see next section). 

The interfacial tension 7' is related to the interfacial 
free energy 7 as [H[ 



7' = 7 + A 



d-, 

dA' 



(11) 



If a liquid is in contact with a dynamic structured wall, 
which can support stress, the interfacial excess free en- 
ergy will vary with the area of the interface. However, in 
this work we consider rigid substrates and structureless 
flat walls, which do not support stress and hence for the 
liquid-wall interface, the interfacial tension will be equal 
to the interfacial free energy validating the use of the PA 
method. For a crystal-wall interface, however, the second 
term in Eq. (|11[) will be a relevant quantity. In this work, 
we will restrict our attention only to the determination 
of the interfacial free energy. 

The interfacial free energy of an inhomogeneous system 
with walls can be defined as a Gibbs excess free energy 
per area, 



7 



^system ^bulk 

A 



(12) 



with G 



system 



NP^AT 



and Gb u ik the Gibbs free energies of the 
inhomogeneous system and the bulk phase of the system, 
respectively. We will use this definition to calculate the 
interfacial free energy using TI. 



B. 7 from PA 

Determination of the interfacial free energy by the PA 
method is only valid if the interfacial tension equals the 
interfacial free energy. This holds, e.g., for interfaces be- 
tween a liquid and a flat wall or rigid substrate. Hence, 
we will use the PA technique to obtain the wall-liquid in- 
terfacial free energy, and compare it with results obtained 
from TI. 

To obtain interfacial free energies from the mechanical 
approach, the local tangential and normal pressure ten- 
sor components have to be computed. There is no unique 
microscopic definition for these local pressure tensor com- 
ponents and different expressions lead to the same value 
for the interfacial tension [24| . Mechanical stability, how- 
ever, requires that the normal component of the pressure 
tensor is independent of the distance from the wall and 
furthermore the two tangential components along the x 
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and y directions are equal to each other. In the literature, 
it is only the Irving and Kirkwood (IK) definition of the 
pressure tensor that satisfies these properties [HI, HE ■ 
According to the IK definition, contributions to the nor- 
mal and tangential components of the pressure tensor 
from any two particles i and j at Zi and Zj , respectively, 
can be written as 



= p(z)k B T 



\i<j 
N p N- 



Z Zi ^ Q ( ?i- z 



\i=i j=i % j 



Z Zi ^ Q ( Z 



(13) 



and 



P^(z) =p(z)k B T 



x6 



2A 



e 



^""^ X ij "t" Vij u pp{ r ) 

z — , 



where is the Heavyside step function, Zi 
p{z) is the local density given by 



N(z) 



{A/2) x Az 



(14) 



and 



(15) 



Here, Az is the bin width used to obtain the pressure 
profiles and N(z) is the number of liquid particles in the 
bin between z and z + Az. This contribution to the local 
pressure tensor is added to all bins between Zi and Zj . It 
is to be noted that the liquid-structured wall interaction 
has no contribution to the tangential component of the 
pressure tensor due to the periodicity of our system in 
the lateral direction [HI [13] ■ 

The contribution to the pressure tensor from the struc- 
tureless walls can also be taken into account by the IK 
method [24|, [25j by considering the walls at Zb and z t 
to be particles of infinite mass. From Eq. flT3"l) . we thus 
obtain 



v i=l 
N 



(16) 



T\£ F lw (z t - Zi)@(z - Zi) \ , 



,4 



with Ff w (z) = —dllf w (z)/dz. 

From Eqs. (fl~3| and ([T4l) . it is clear that if two particles 
in a bin are located on the same side of z, their contri- 
bution to the local pressure tensor cannot be taken into 
account by the IK method. To minimize the number of 
such cases, we must choose the bin width to be compara- 
ble to the shortest distance between between the particles 
in the z direction. On the other hand if the bin width is 



too small, there will be larger fluctuations in the pressure 
tensor and the average must be taken over many more 
configurations to get a smooth profile, thus increasing 
the computational time. In our simulations we choose a 
bin-width of Az = 0.05. 

Equation (JTUJ) being the difference between two sim- 
ilar numerical values is subject to large relative errors. 
Moreover, at large densities near the wall, the density 
and pressure profiles show rapid oscillations and hence 
resolving them with high precision requires a huge com- 
putational effort. Below, the accuracy of the PA method 
is studied in detail via a direct comparison to the data 
obtained from TI. 



C. 7 from TI 

In a TI, the free energy of a state of interest is com- 
puted with respect to a reference state 2l[ . A parameter 
A, which couples to the interaction potential, 



is gradu- 
ally changed such that the reference state is reversibly 
transformed into the final state of interest. 

To calculate the interfacial free energy of the LJ sys- 
tem in contact with a structured wall, the TI scheme is 
carried out in two steps. In the first step, a bulk LJ sys- 
tem without walls and periodic boundary conditions in 
all directions is reversibly transformed into a LJ system 
in contact with a structureless flat wall along the z di- 
rection. In the second step, the flat wall interacting with 
the LJ system is reversibly transformed into a structured 
wall. To ensure reversibility of the thermodynamic path, 
periodic boundary conditions are applied in x, y and z 
direction. Calculating the free energy change in the two 
steps yields the required interfacial free energy. 

To obtain 7 for a hard-sphere system in contact with 
a hard structureless wall via TI, Heni and Lowen [29[ 
have used a scheme, where a bulk hard sphere system 
is reversibly transformed into a system interacting with 
an impenetrable hard wall. In this work, we generalize 
the scheme of Heni and Lowen to continuous wall po- 
tentials. To this end, the wall potential is parametrized 
by a parameter A such that the wall changes smoothly 
from a penetrable to an impenetrable wall as A increases. 
The following parametrization of the wall potential is 
adopted: 



Wf w (A, z) =A 2 4e v 



z + (1 - A)z cw 

> 6 

'pp 



z+(l- X)z c 



(17) 



x w(z). 



Figure [U shows the parametrized wall potential at dif- 
ferent values of A. At A = a bulk LJ system can freely 
cross the boundaries. For small values of A the barrier 
height at z — is of the same order as k B T and the LJ 
particles can penetrate the barrier. As A increases, the 
wall becomes more and more impenetrable and finally an 
impenetrable WCA wall is obtained at A = 1. 
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Since, the interfacial excess free energy of the LJ sys- 
tem in presence of walls is calculated with respect to a 
bulk LJ crystal, it is important that the bulk density is 
maintained as the parameter A is varied during the trans- 
formation. This is particularly important for a LJ liquid 
close to coexistence, since an increase in the bulk density 
in the presence of walls could lead to a precrystallization 
of the bulk liquid during the transformation, thus mak- 
ing it irreversible. A constant bulk density also ensures 
that our system is large enough such that there are no 
mutual influences between the walls on either side of the 
bulk L J system. To maintain a constant bulk density one 
must keep the normal pressure Pn constant and change 
the volume, as is achieved by carrying out simulation in 
the NP N AT ensemble. 

The system Hamiltonian now depends on A and is 
given by 



A is 



N p N p N p 



i=l 



i—1 



N„ 



(18) 



2ju fw (A, z = Zi-Z) 



and thus the partition function can be written as 

H(r, p, A) + P Z AL Z 



0(A) 



1 

h 3N M 



exp 



xAdL z dp N dp N 



(19) 



The derivative of the Gibbs free energy with respect to 




FIG. 1: (Color online) Variation of the WCA wall potential as 
a function of A during the transformation of the bulk LJ liquid 
or crystal into an impenetrable flat wall interacting with the 
LJ system. From no wall at A = 0, we have a wall with a 
finite barrier at small values of A. With increasing A, the wall 
becomes more and more impenetrable. At A = 1, there is an 
impenetrable wall represented by the WCA potential. 



dG(X) 

dX 



0(A) 



dQ{\) 



dX 



(20) 



where the angular brackets denote the ensemble average 
at a particular value of A in the NPrAT ensemble. 

The Gibbs free energy difference between the two ini- 
tial and final state can then be obtained as 



AG = G(X = 1) - G(X = 0) 
- 1 /dH{X)\ 



f 


dG(xy 




<9A 



dX . 



dX (21) 



(22) 



To compute AG from molecular simulations, indepen- 
dent simulations runs are carried out at N\ discrete in- 
tervals between A = and A = 1. Alternatively, one 
can also calculate the free energy difference in a single 
simulation by varying A step by step such that the final 
configuration at a value of A = A^ is the initial configura- 
tion for the next value at A = Aj+i. In both methods, the 
system is equilibrated at each A = A<j, and then the time 
average of the quantity dH(X)/dX is calculated. The nu- 
merical integration of Eq. (|2"2"|) is carried out using the 
trapezoidal rule: 

AG = E 2 ^ dH ' dx ^ + ( dH / dx ^+i] ( A *+i - A *) ■ 

i=l 

(23) 

The partial derivative of H(X) with respect to A is given 
by 



dH(X) du iw (z, A) 



dX 



dX 



X 



z+(l - X)z c 



Ufw(A, z)-\ 
-Ufw(z, A) 



(24) 



The above TI scheme leads to a wall which is not fully 
impenetrable and hence does not correspond to the de- 
sired state of interest at the end of the integration path. 
While the LJ particles cannot cross the wall at A = 1, 
two particles near the boundary but on opposite sides of 
the wall can still interact with each other. To overcome 
this problem another TI step is carried out to bring the 
system to a state where the LJ particles are in contact 
with a fully impenetrable wall excluding such spurious 
interactions. This is achieved by parametrizing u pp (r) 
by a factor fj,, 



<<pp 



= u pp( r ij) + (l- V)u 



(n 3 ) (25) 



where, upp vij) denotes interaction between LJ particles 

(2) 

on same side of the wall, while C/ pp (ry) corresponds to 
interaction between particles near the boundary but on 
opposite sides of the wall, i.e. the separation between par- 
ticles is greater than L z /2. At /i = 0, all such spurious 
interactions are taken into account. As /i increases such 
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interactions are reduced by the factor 1 — /U and finally 
at fi = 1, these spurious interactions are completely ne- 
glected. The fi dependent Hamiltonian for this step can 
be written as, 



i—l 

E E [^w + a-M^M 



(26) 



2—1 



^ Wfw(A = 1, z = Zi - Z) . 

i=i 

The thermodynamic integrand in this step is 

dH/dfi = du pp (r,n)/d(i = -«^(r«). (27) 
and thus the free energy difference can be expressed as 



AG 



fw—>fw* 



(du pp (r,fx)/dfj,)^dp,. 



(28) 



Our simulations showed that the contribution of 
AGf w _j.f w * is very minor, i.e. about 0.1% of AG from 
Eq. (j2"3")l and hence can be neglected. 

Using Eqs. (JT2J), ([2"T]l. and pi)), the interfacial free en- 
ergy of a LJ system interacting with a flat wall can be 
written as 



G 



fw 



Gbuii 



AG 



bulk— >iw 



A 



with 



AG bu 



lk— >fw 



du[ w (z, X) 
dX 



dX. 



(29) 



(30) 



In the second step of our TI scheme, the flat wall is 
reversibly transformed into a structured wall in contact 
with the LJ system. During this change, the flat walls 
are positioned at the same location as the structured wall 
layer closest to the LJ liquid or crystal and there is no 
interaction between the flat and structured walls. The 
transformation from flat walls to structured walls is ac- 
complished by parametrizing the wall potential as: 



U wa xi(nj,X) = (1-A) 2 y^u fw (2 



-Z)+A 2 ^u pw (rij) . 

(31) 



Now, the A-dependent Hamiltonian is 

2m 

i=i 



^p- A )-E^7-Pi +E E u pp( r «)+ 

1=1 j=i+i 

(1 - A 2 ) ^ u fw( z = z i - z - A) + A 2 ^2 E M P w ( ri j) 

i—l i—l j—1 

(32) 



and the derivative of the Hamiltonian with respect to A 

dH _ dt/ wa n 
dX ~ dX 

=2 



(a - i)y\ fw (z = Zi—z) + Ay^y^M PW (rij) 

i i=l j = l 

(33) 



So, finally the interfacial free energy of the LJ system in 
contact with a structured wall (sw) is given by 

G S w — Gbulk AGbulk— >fw , AGf w _j. sw / nA \ 

7wc = ; = j + - A (34) 



A 



with 



AG- 



fw- 



V ^wall(A) 

dX 



A 



dX. 



(35) 



D. Simulations 

To integrate the equations of motion, the velocity 
form of the Verlet algorithm was used with a time step 
r = 0.005 and, to maintain constant normal pressure, 
the Andersen barostat algorithm was chosen. Peri- 
odic boundary conditions are employed in the x, y and z 
directions for the first step of the TI method where flat 
walls are considered. In the second step periodic bound- 
ary conditions arc only used along the x and y direc- 
tions. The PA simulations are carried out with periodic 
boundary conditions only along the x and y directions. 
The temperature was kept constant by drawing every 200 
steps the velocity of the LJ particles from the Maxwell- 
Boltzmann distribution at the desired temperature. 

During the NP^AT simulations, the position of the 
flat or structured walls must be modified keeping the nor- 
mal pressure Pn constant. To ensure this, the flat walls 
are treated as particles of infinite mass and, at each time 
step, the wall position 2f w is rescaled according to 



Z(t + At) = Z(t) x L z (t + At)/L z (t) 



(36) 



Note that this method is similar to the "fluctuating wall" 
method proposed by Lupowski and van Swol [4l| , main- 
taining a constant normal pressure in a MC simulation 
of LJ particles in presence of a structureless wall. 

When a rigid structured wall interacts with the L J sys- 
tem, the wall particles must not change their positions 
relative to each other, thereby changing the wall density. 
To circumvent this problem, the center of mass of the 
wall is changed at every time step according to Eq. (|36p . 
The position of the individual particles of the wall is then 
shifted such that they are at the same relative distance 
from the center of mass as at the beginning of the simu- 
lation. 

To calculate 7 w i we consider systems of 4000 particles. 
The structured walls contain between 200-1200 particles, 
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depending on the orientation and the density of the wall. 
The total surface area of the simulation cell is about 
A = 200, yielding a length along the z direction of about 
L z = 65 at the various wall-liquid interaction strengths 
e pw and structured wall densities p w . 7 w i is computed at 
a normal pressure of Pn — 3 and temperature T = 2. At 
the start of the simulation, the LJ particles were placed 
on ideal fee lattice sites and the walls were inserted si- 
multaneously. Then the system was allowed to melt and 
equilibrate at the desired pressure and temperature, be- 
fore the calculations were performed. 

To test for the presence of any finite size effects, we 
also performed simulations with up to 12000 particles 
and a total surface area of about A = 340, but obtained 
identical results compared to the simulations carried out 
with the smaller system size. This shows that systems of 
4000 particles are large enough to avoid finite size effects 
in the calculation of interfacial free energies. 

From previous works pertaining to hard sphere sys- 
tems, it is well known that the (111) orientation of the 
crystal in contact with a planar hard wall (or a soft WCA 
wall) gives the lowest interfacial tension as compared to 
the (100) or (110) orientations [42j |. At small undercool- 
ings, the hard sphere fluid freezes into the (111) crystal 
near the wall [43|]. Hence, we obtain interfacial free ener- 
gies only for the (111) orientation of the fee crystal phase 
in contact with the walls. Unlike the liquid, the crystal 
has a long-range order and, in order to prevent deforma- 
tion of the crystal, the system size must be compatible 
with this order. 

For the determination of 7 WC , systems of 7056 particles 
and area around A = 450 are considered. The number of 
structured wall particles ranges from 800 to about 1200, 
depending on the different wall densities. Only the (111) 
orientation of the crystal in contact with the (111) orien- 
tation of the structured wall along the interface normal 
was considered. The corresponding simulations to obtain 
the interfacial free energy of a crystal in contact with a 
flat wall are carried out with 3960 particles with an area 
of around A = 200. Simulations were also carried out 
with a system size of 6006 particles and a total area of 
A = 300 and there was only a marginal deviation (< 1%) 
in the value of 7 WC as compared to the smaller system. 

For comparing results obtained by our approach with 
that of Tang-Harris [l6| and Grzelak-Errington [36|, we 
performed simulations for their system with 4000 liquid 
particles and 392 structured wall particles at the tem- 
perature T — 0.9. We considered a lateral system size of 
10 x 10 and the length of the box along the z direction was 
kept at 60.5134 to obtain a bulk liquid density of 0.661, 
the value reported by Tang and Harris [l6| for their simu- 
lations. Simulations were performed at this fixed density 
in the NVT ensemble. With this system size, the finite 
size effects were negligible. The liquid in contact with 
the flat wall [Eq. (gj)] was used as the reference state to 
calculate 7 w i for the liquid in contact with the structured 
wall at the same bulk density and temperature. 

A simulation at constant normal pressure leads to fluc- 



tuations of the length of the simulation cell in the z di- 
rection, L z . However, in order to compute the density 
and pressure profiles necessary for the PA method, it is 
more suitable to keep L z constant. Hence, to obtain 7 w i 
via the PA method, we first equilibrate the system in the 
NP^AT ensemble for 5 x 10 5 time steps. After equilib- 
rium is reached, the simulations continue for 4.5 x 10 6 
time steps, from which the average length of the box in 
the z direction is calculated. L z is set to this average 
value and the particle coordinates are rescaled by the 
factor (L z ) / L z (t{) , U denoting the time at the end of this 
equilibration run. An equilibration run is then carried 
out in the NVT ensemble for 5 x 10 5 time steps and the 
final production run consists of 4.5 x 10 6 steps, when we 
accumulate data for the density, energy, temperature and 
pressure profiles every 5 time steps, averaging the profiles 
over 9x 10 5 sample configurations. In our simulations, we 
observe a drift of 0.5 — 2.5% in the normal pressure profile 
from the given external pressure Pn. This drift can be 
reduced by averaging the length of the box for a longer 
simulation time or over a large number of realizations. 

To calculate the interfacial free energy via TI, we used 
around 40 intervals between A = and A = 1 to nu- 
merically compute Eq. (1231) . Independent equilibration 
runs were carried out at each value of A, in the NP^AT 
ensemble for about 5 x 10 5 — 1 x 10 6 time steps. After 
the completion of the equilibration run, production runs 
were performed for 5 x 10 5 steps in order to accumulate 
data. The same TI scheme and simulation procedure has 
been adopted to determine the interfacial free energy of 
the system investigated by Tang-Harris [l6[ and Grzelak- 
Errington (3(| , but in the NVT ensemble at a fixed liquid 
density. 




10 20 30 40 50 60 
z 



FIG. 2: (Color online) Density profile of the system configu- 
ration for different values of A at the temperature T = 2.0, 
the normal pressure Pn = 3.0, and the wall-liquid interaction 
strength e w = 1. The corresponding density profile for a liq- 
uid in contact with the (100) orientation of a structured wall 
of density p w = 1.371 at e pw = 1.0 is also shown. The inset 
shows the density profiles in the bulk region on a magnified 
scale. 
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FIG. 3: (Color online) (a) {duf vr (z, X)/dX) as a function of A computed from simulations at Pn = 3 and T — 2 for the liquid and 
at Pn = 3.0 and T = 0.5 in case of the crystal. To determine -y w i, 4000 liquid particles were enclosed in a simulation box of area 
A — 200 with the wall-liquid interaction strength e w = 1. Corresponding simulations to calculate 7 WC were carried out with 
3960 particles and and area A = 235.22. (b) {d[U-mjx\/dX), corresponding to transformation of the flat wall into a structured 
wall. In case of the liquid, the structured wall consists of 392 particles rigidly fixed to fee lattice sites, with the (100) orientation 
of the wall facing the liquid. For the crystal, the structured wall consisted of 432 wall particles with the (111) orientation of the 
wall in contact with the crystal. The density of the structured wall was p w = 1.371 for the liquid (p rm = 0.647 for the crystal) 
and the wall-particle interaction strength for the liquid -wall simulations was kept at e pw = 1, while for the crystal e pw = 0.5. 
Other parameters are same as in (a). 



For our TI method to be valid, there must be a bulk 
region unaffected by the wall. Figure [5] shows the density 
profile of liquid in contact with the parametrized flat wall 
represented by Eq. (|17|l at various values of A, with the 
wall-liquid interaction strength e w = 1. Also shown in 
Fig. [2] is the density profile of the liquid in contact with 
the (100) orientation of a structured wall of density p w — 
1.371 and with interaction strength e pw = 1. The inset 
shows a magnified view of the density profiles in the bulk 
region. Clearly, all the density profiles overlap with each 
other indicating that the bulk region is unaffected by the 
walls. 



Figure shows the thermodynamic integrand 

(du[ w (z, A) / d\) as a function of A, during the transforma- 
tion of a bulk liquid (crystal) to a confined liquid (crys- 
tal), interacting with flat walls. The integrand is smooth, 
thus allowing for an accurate determination of the inter- 
facial free energy. Figure \3jp shows the integrand as a 
function of A for the second step of the thermodynamic 
integration when the flat wall is transformed into a struc- 
tured wall. The integrand is always negative, implying 
that the interfacial free energy of a LJ liquid (crystal) in 
contact with a rigid structured wall is smaller than for 
the case where the liquid (crystal) is in contact with a 
structureless flat wall. 



IV. RESULTS 

A. 7 w i 

Using TI, we first determine the liquid-flat wall interfa- 
cial free energy 7 w i at several temperatures and pressures. 
In Fig. [H we plot the data obtained from TI along with 
the estimate for 7 w j from the PA technique, as a func- 
tion of temperature and pressure, respectively. The error 
bars in 7 w i obtained from TI are smaller than the symbol 
size and hence are not reported. It is evident from Fig. 0] 
that there is good agreement between the two methods 
within the statistical error. However, Fig. 0^ shows that 
7wi obtained from TI is smoothly varying, while the PA 
data is less systematic. Relative differences between the 
two methods are between 0.3 and 1.8%. 

This small disagreement between the two methods is 
due to the large fluctuations in the local pressure pro- 
files as obtained from Eqs. f| 13[) and (TT4")) (see Fig. [5]). 
The inset in Fig. [5] clearly shows the large fluctuations 
in the normal pressure profile near the wall and in both 
the normal and tangential pressure profiles in the bulk 
region. Since Eq. (1101) represents the difference between 
two pressure profiles: Any lack of precision in the numer- 
ical measurements magnifies the relative error. The TI 
data is more accurate and less computationally expen- 
sive compared to what would be required to obtain more 
precise values from the PA method. 

The liquid-flat wall system can now be used as the ref- 
erence system to calculate the interfacial free energy of 
the liquid in contact with a rigid structured wall. Some 
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FIG. 4: (Color online) (a) Interfacial free energy of liquid in 
contact with the flat wall, 7„i, as a function of temperature. 
NPnAT simulations were carried out at a normal external 
pressure Pn = 3.0, with 4000 liquid particles and total surface 
area of A = 200. The wall-liquid interaction strength is e w = 
1. Filled squares correspond to data obtained from TI, while 
estimates from PA are represented by open circles with error 
bars. Uncertainty in data computed by TI is less than the 
symbol size, (b) 7 w i as a function of pressure at T = 2.5. 
Other parameters and symbols representing the TI and PA 
data are same as in (a). 



10 
9 
8 
7 
6 
5 
4 
3 
2 
1 




- P N (z) 3.08- 
P T (z) 

£3.04- 

3- 
L 

J 


... 


) 20 40 60 

| 





10 20 30 40 50 60 
z 



FIG. 5: (Color online) Normal and tangential components of 
the pressure profiles of liquid in contact with a flat wall at 
P N = 3 and T = 2, with 4000 liquid particles. The wall liquid 
interaction strength, e w = 1. Inset shows the pressure profiles 
on a magnified scale close to the magnitude of the external 
pressure Pn. 
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FIG. 6: (Color online) a) A two-dimensional projection onto 
the zy plane of a sample configuration of liquid in contact 
with the (100) orientation of the structured wall at Pn = 3 
and T — 2. 4000 liquid particles are enclosed in a simulation 
box of area A = 200. The density of the structured wall is 
p w = 1.371. The wall-liquid interaction strength is set to 
e P w = 1. b) Density profile of the liquid, averaged over many 
configurations, showing pronounced layering at the structured 
wall-liquid interface. 



properties of the structured wall such as the wall-liquid 
interaction strength, density of the structured wall and 
its orientation along the interface will affect the inter- 
facial free energy and consequently the wetting behav- 
ior of the liquid. We will investigate 7 w i for the (100), 
(110) and, (111) orientations of the structured wall in 
contact with the liquid at different wall-liquid interac- 
tion strengths. Effects of density of the structured wall 
on the interfacial free energy will also be studied. Un- 
less otherwise indicated, the external pressure is set to 
Pn = 3 and the temperature to T = 2. 

Figure [6] shows a sample configuration of the liquid in 
contact with the (100) orientation of a structured wall 
in the zy plane and the corresponding density profile. 
We observe layering of the particles near the wall. Away 
from the walls a bulk region forms, where the density is 
constant. 

In Fig. [7^, 7wi is displayed as a function of e pw , as ob- 
tained from the PA and TI method. The error bars in the 
PA method were calculated from 2 — 3 realizations. For 
TI, the error bars are smaller than the size of the sym- 
bols and hence they are not reported. We observe that 
7wi decreases with e pw , which is in agreement with previ- 
ous studies carried out using the mechanical route [HI, EH 
and other thermodynamic methods p^ . [36j . High values 
of e pw represent stronger attraction between the wall and 
liquid particles. This reduces the free energy needed to 
move the liquid from the bulk to the surface resulting in 
a lower interfacial free energy. Data from the two meth- 
ods agree qualitatively but the percentage difference of 
the PA results with respect to the TI data increases at 
higher values of e pw . 

Apart from the fluctuations in the local pressure pro- 
files, the strong layering near the interface at large e pw 
also reduces the numerical accuracy of the PA method. 
Figure [TJs shows the liquid density profiles at various in- 
teraction strengths along with the corresponding profile 
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FIG. 7: (Color online) a) Interfacial free energy of liquid in contact with the (100) orientation of the structured wall vs. e pw . 
All parameters are the same as in Fig. [6] Open squares with solid line represent data from the TI method and open circles with 
dashed line represent estimates yielded by the PA method. Uncertainties in data from TI are less than the symbol size. The 
wall-liquid interaction strength is e pw = 1. b) Liquid density profiles at different wall-liquid interaction strengths. The inset 
shows the profiles corresponding to the (100), (110) and (111) orientations of the structured wall in contact with the liquid at 
e P w = 1. For comparison the density profile for liquid in contact with the flat wall is also shown in both a) and b). All system 
parameters are same as in a). 



in presence of a flat wall. The first peak in the density 
profile corresponding to the flat wall occurs at a greater 
distance from the wall as compared to the peaks arising 
out of the liquid-structured wall interaction. This can be 
attributed to the purely repulsive flat wall, which pushes 
the liquid further away from the walls compared to the 
structured walls. 

The interfacial free energy of the crystal-melt interface 
is influenced by the orientation of the crystal in contact 
with the melt [l|- Similarly, it might be expected that 
different orientations of the structured wall in contact 
with the liquid will affect the wall-liquid interfacial free 
energy. In the inset of Fig.[7]D, we plot the density profiles 
near the wall for the (111), (110) and (100) orientations of 
the structured wall in contact with the liquid at e pw = 1 . 
The density of the structured wall p w = 1.371 and, the 
lateral dimensions of the system corresponding to the 
(100), (110) and (111) orientations of the wall are 10x10, 
10 x 10.102 and 9.623 x 10.102, with 392, 420 and 330 
wall particles, respectively. Figure [7p shows the layering 
of the density profile to be most pronounced for the (111) 
orientation owing to the closely packed atoms exerting 
a greater repulsive force on the liquid. In contrast, the 
layering for the (110) orientation is much less pronounced 
and the first peak in the density profile also occurs closer 
to the wall as compared to the (111) or (100) orientations. 

In Table I, we report y w i, obtained from the PA and 
TI methods, at various e pw , for the (100), (110) and 
(111) orientations of the structured wall in contact with 
the liquid. In general, we find 7 w i corresponding to the 
(111) and (100) orientations of the wall to be larger com- 
pared to the (110) orientation. This can be attributed 
to the stronger repulsive forces exerted on the liquid by 



the more close packed (111) and (100) planes, as com- 
pared to the more loosely packed (110) plane. Also, rela- 
tively better agreement is observed between the PA and 
TI methods for the (100) and (111) orientations of the 
wall as compared to the (110) orientation. 

Simulations were also carried out for large system sizes 
at e pw = 1, with up to 12000 particles, and large sur- 
face area and wall separations. However, no systematic 
change was observed as compared to the smaller system 
size. Clearly, to obtain accurate values of the interfacial 
free energy, the pressure profiles need to be determined 
with far greater numerical accuracy. This has also been 
recently pointed out by D. Deb et al. [H[ for hard sphere 
systems. Computing the pressure profiles with high pre- 
cision is computationally expensive and since accurate 
values can be obtained by the TI method with much less 
computational effort, use of the PA technique seems to 
be unjustified. In the remainder of the discussion on our 
model, we report results obtained with the TI method 
only. 

The density of the structured wall will also have an 
impact on the wetting behavior of the liquid in contact 
with it. We have carried out simulations at several den- 
sities p w corresponding to different lattice constants of 
the ideal fee lattice structure of the wall. In Fig. [5£i, 
we report the TI results for 7 w i as a function of p w at 
three different e pw 's. At e pw = 1, the interfacial free en- 
ergy decreases with the density of the wall. The larger 
number of wall particles at greater densities exert strong 
attractive forces on the liquid, reducing the interfacial 
free energy. At extremely large densities the interfacial 
free energy becomes negative indicating that the liquid 
completely wets the wall. No data for very low densities 
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0.10 


2.115 


2.165±0.034 


1.865 


1.902±0.025 


2.149 


2.189±0.010 


0.25 


2.039 


2.067±0.034 


1.808 


1.860±0.033 


2.056 


2.073±0.004 


0.50 


1.742 


1.775±0.012 


1.521 


1.587±0.003 


1.740 


1.806±0.002 


0.75 


1.364 


1.343±0.007 


1.144 


1.252±0.003 


1.348 


1.335±0.002 


1.00 


0.929 


0.862±0.012 


0.705 


0.868±0.012 


0.906 


0.902±0.012 



TABLE I: Interfacial free energy 7 w i at different wall-liquid interaction strengths, for the (100), (110) and (111) orientations 
of the structured wall in contact with the liquid. Data computed from both TI and PA are shown. Simulations are carried out 
at Pn = 3 and T — 2. The density of the structured wall p w = 1.371. 
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FIG. 8: (Color online) a) Interfacial free energy of a liquid in contact with a structured wall vs. density of the wall at several 
wall-liquid interaction strengths e pw . Results were obtained from TI. Simulations were carried out at a constant normal pressure 
Pn = 3 and temperature T = 2.0. b) Density profiles of the liquid in contact with the structured wall at various densities p w . 
The wall- liquid interaction strength is e pw = 1. Other parameters are same as in a). 
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FIG. 9: (Color online) 7„i vs. e pw for the model specified by 
Eq. ([6]). Squares correspond to our results computed via TI. 
Filled circles represent data from the studies of Tang-Harris 
[l6T |. while diamonds and filled right triangles are the data 
obtained by Grzelak-Errington [36j using the one-wall and 
two-wall approaches, respectively. 



of the structured wall are shown in Fig. [5^ since the liq- 
uid particles penetrate the wall at low densities and the 
interfacial region is no longer well defined. 

At e pw = 0.5, 7 w i shows a weak maximum and at large 
wall densities decreases with p w more gradually as com- 
pared to the situation when e pw = 1. At still lower 
wall-liquid interaction strength (e pw = 0.1), 7 w i has a 
weak dependence on the density of the structured wall 
and remains almost constant in the range of p w shown in 
Fig. [8k, 

In Fig. [SJd, we show the density profiles of the liquid 
corresponding to several densities of the structured wall 
p w at e pw = 1. It is observed that the layering gets 
more pronounced and the first peak in the density profiles 
occurs further away from the walls as p w increases. A 
similar behavior of the density profile was observed when 
increasing e pw at fixed p w . The variation in j w \ and the 
nature of the density profiles indicate that increasing p w 
at e pw = 1 has a similar effect on 7 w i as increasing e pw 
at a fixed value of the wall density. 

Finally, to compare results from our TI technique with 
those obtained by other methods, we consider the model 
defined by Eq. (J5|), which was first studied by Tang and 
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Harris using a PA technique [l6j]. Later Grzelak and 
Errington [36| utilized GCTMMC simulations to obtain 
free energy profiles of the same system over a wide range 
of densities, with the fluid confined by a structured wall 
on one side and a hard wall at the other side or the fluid 
confined between two identical structured walls. Both of 
their approaches with one or two structured walls lead 
to the same results within the statistical errors. Using 
Tl in the NVT ensemble, we computed the interfacial 
free energy of the same system, with the liquid confined 
by two identical structured walls. In Fig. |H1 our re- 
sults are reported along with data from the two previous 
works [IH, [36| . Data obtained by Tang-Harris systemat- 
ically deviate from our estimates of 7wi &s 6pw increases. 
Their data also has a large statistical error. However, our 
predictions are in good agreement with those of Grzelak 
and Errington. 
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FIG. 10: (Color online) Crystal-flat wall interfacial free en- 
ergy 7 w i (diamonds) and liquid-flat wall interfacial excess free 
energy 7 WC (squares) as a function of temperature. NP^AT 
simulations were carried out at Pn = 3 and the wall-liquid 
interaction strength e w = 1. The (111) orientation of the 
crystal was considered in determining 7 WC . Inset corresponds 
to 7 w i and 7 WC as a function of the interaction strength e w at 
coexistence: Pn = 3.0 and T = 0.8. Symbols are same as in 
the main graph. 



We will compute the crystal-wall interfacial energy by 
the TI method only, since the crystal can support stress 
and hence the interfacial tension and interfacial free en- 
ergy are not the same, thus invalidating the applicabil- 
ity of the PA method [23j. In performing simulations of 
crystal in contact with walls on both sides, the number 
of particles N must be chosen such that it is compatible 
with the long range order of the crystal. This is in con- 
trast to a liquid, where choosing a large enough N yields 
a sufficiently large system and the two walls on either 
side of the liquid do not influence each other. However, 



the crystal has a long-range order and merely choosing a 
large N may not necessarily be commensurate with this 
order. Such an incommensurate N is associated with 
long range elastic distortion, that propagates from one 
wall to the other leading to an inaccurate value for 7 WC . 
This had already been pointed out by D. Deb et al. [32| 
who studied the interfacial free energy of a hard-sphere 
crystal confined between softly repulsive walls described 
by the WCA potential. 

As specified earlier, we restrict our attention to the 
close-packed (111) orientation of the crystal in contact 
with a flat wall along the z axis. To evaluate 7 WC , a 
bulk fee crystal with the (111) orientation along z axis 
is simulated in the NPT ensemble. Periodic boundary 
conditions are employed in all directions to determine 
the average equilibrium lattice constant and hence the 
density of the crystal. A fee crystal with this density 
was chosen as the initial configuration for the TI sim- 
ulations to compute 7 WC . The length of the simulation 
box was chosen such that an integer number of unit cells 
along the x, y and z directions adapted exactly into the 
simulation box. Then, independent simulations in the 
NP^AT ensemble were carried out at each value of the 
A parameter during the two-step TI scheme. In Fig. [TOl 
we plot 7 WC for crystal in contact with a flat wall, as a 
function of temperature up to the coexistence tempera- 
ture at Pn = 3.0. For comparison, 7 w i is also plotted 
at the same pressure. Similar to 7 w i, 7 WC decreases as a 
function of temperature. 

To predict the wetting behavior of the crystal in con- 
tact with the wall at crystal-liquid coexistence, one needs 
7d in addition to 7 w i and 7 WC . However, without knowl- 
edge of 7d, it is still possible to predict whether the crys- 
tal will completely wet the wall (8 C = 0°) or do so only 
partially (0 < 8 C < 180°). To this end, simulations were 
carried out at coexistence (Pn = 3.0, T = 0.8) for the 
bulk liquid and crystal in contact with a flat wall at var- 
ious interaction strengths between the wall and the bulk 
liquid or crystal. The data is reported in the inset of 
Fig. [TO] No wetting layer was observed near the walls 
during wall-liquid simulations, allowing for a determina- 
tion of the wall-liquid interfacial free energy directly at 
coexistence. We find that 7 WC > 7 w i, showing that there 
is incomplete wetting of a fiat wall by the (111) orienta- 
tion of the LJ crystal and that the contact angle can be 
varied by changing e pw . 

While 7d has not been determined in this work, David- 
chack and Laird [3 Ef| , obtained the crystal- liquid in- 
terfacial free energy at coexistence for a similar LJ model. 
The bulk liquid and crystal densities at the coexistence 
temperature T = 0.809 for their model is same as for our 
system at P N = 3.0 and T = 0.8. At T = 0.809, they 
obtained 7 c i = 0.428 ± 0.004. Since the LJ model used 
in this work is not very different from their potential, it 
is safe to assume that the crystal-liquid interfacial free 
energies will not be far apart. Using j c \ = 0.428 ± 0.004 
and the values of j w \ and 7 WC used in this work, we ob- 
tain contact angles of 97.4°, 103.8° and 113.6° respec- 
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tively, signifying partial wetting of the flat wall. This 
is in contrast to the hard sphere case, where the (111) 
orientation of the crystal led to complete wetting of the 
wall [32[ . Such a situation of incomplete wetting will fa- 
cilitate the study of heterogeneous nucleation of a crystal 
droplet at a wall-liquid boundary, and enable us to test 
the predictions of classical nucleation theory. 
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FIG. 11: (Color online) Crystal-structured wall interfacial free 
energy, 7 WC , as a function of the structured wall density, for 
the (111) orientation of the crystal in contact with the (111) 
orientation of the structured wall. NPnAT simulations were 
carried out at a normal pressure and temperature Pn = 3 
and T = 0.5 respectively, with the crystal-wall interaction 
strength e wc = 0.5. The inset shows 7 WC as a function of e pw 
for a structured wall density p w = 0.647, other parameters 
remaining the same. 

Having obtained the flat wall-crystal interfacial free 
energy, we can now compute the structured wall-crystal 
interfacial excess free energy 7 WC . We choose to investi- 
gate the (111) orientation of the crystal in contact with 
the (111) orientation of the structured wall. To obtain 
7 WC , commensurate surfaces of the wall must be in con- 
tact with the crystal on both sides. We know that in the 
fee structure, there is an ABC ABC ABC... stacking of 
the lattice planes along the (111) orientation. The same 
order of the planes must be kept for the crystal plane 
in contact with the wall. For example, the the following 
stacking of the planes, 



A W B W C W A W B W C W A C B C C C A C B C C C ... 
A C B C C C A W B W C W A W B W C V/ . 



(37) 



is commensurate. However, a stacking of the planes in 
an incommensurate manner such as 



A W B W C W A W B W C W A C B C C C A C B C C C ... 
A c B C C C C W B W A W C W B W A W 

will lead to long range deformation of the crystal. 



(38) 



In Fig. [TT] and its inset, we plot 7 WC as a function of 
the structured wall density and, in the inset, as a func- 
tion of the wall-crystal interaction strength e pw . Similar 
to the liquid case we find that the interfacial free en- 
ergy decreases with e pw due to the stronger attraction 
between the crystal and the wall. Unlike the liquid case, 
Fig. [TTJ shows that while the main trend for 7 WC is to 
increase with decreasing density of the structured wall, 
there is a sharp dip when the density of the wall equals 
the density of the crystal. This is easy to understand, 
since less energy will be needed to create an interface, 
when the structured wall has the same structure as the 
crystal than when there is a mismatch between the wall 
and crystal structures leading to a relatively unfavorable 
interaction between them. 



V. CONCLUSION 

We propose a thermodynamic integration (TI) scheme 
to compute interfacial free energies of liquids or crystals 
in contact with flat or structured walls from molecular 
dynamics simulation. In this work, this scheme has been 
applied to Lennard- Jones systems, but it can be easily 
generalized to other interaction models. The implemen- 
tation of our method is simple, and, as demonstrated 
above, our method provides reliable and accurate esti- 
mates of 7 w i and 7 WC that enter in Young's equation JJJ. 
In particular for structured walls (substrates) , to the best 
of our knowledge, there are no simulation studies cal- 
culating the substrate-crystal interfacial free energy 7 WC . 
Most of the previous simulation works on structured walls 
[IB - fill . I2R \2l\ have been limited to the calculation of the 
interfacial free energy 7 w i using the integration over the 
pressure anisotropy (PA). The PA method, however, does 
not give reliable results in general, and, in contrast to 
our TI scheme, it is not applicable to substrates that can 
support stress (such as structured walls where the wall 
particles are allowed to move and are thus not fixed to 
their ideal lattice positions, see discussion above). There- 
fore, the TI scheme proposed in this work can be consid- 
ered as a novel approach to obtain accurate values for 
substrate-liquid or substrate-crystal interfacial free ener- 
gies and thus it will be useful in studies of wetting and 
nucleation problems. 
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